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1.   INTRODUCTION 

In  this  paper  we  will  introduce  a  simple  method  for  obtaining 
a  stationary  sequence  of  dependent  random  variables  having  a  specified 
marginal  distribution  and  correlation  structure  (second-order  joint 
moments).   One  advantage  of  the  model  is  that  the  specification  of  these 
two  aspects  of  the  model  is  independent.  Another  advantage  is  that  the 
sequence  is  obtained  as  a  very  simple  transformation  of  a  sequence  of 
independent  random  variables.  The  model  is  analogous  to  models  for 
dependent  sequences  of  exponential  random  variables  introduced  in  Jacobs 
and  Lewis  (1977)  and  Lawrance  and  Lewis  (1977). 

We  will  now  define  some  quantities  which  will  be  used  throughout 
the  paper.  Let  {Y  }  be  a  sequence  of  independent  random  variables  taking 
values  in  a  discrete  space  E  each  having  the  distribution  it.     Let 
{U  )   and  {V  }  be  independent  sequences  of  independent  random  variables 
taking  the  values  0  and  1  with 


(1.1)  P{Un  =  1}  =  p    and    P{Vn  =  1}  =  p 


for  fixed  0  <  0  <  1  and  0  <  p  <  1.  Let  {Sn}  be  a  sequence  of 
independent  identically  distributed  random  variables  taking  the  values 
0,1.2,  ...,N  with  distribution  F,  where  N  is  a  fixed  integer. 

The  most  general  case  which  we  will  consider  here  is  a  sequence 
of  random  variables  {X  }  which  is  formed  according  to  the  probabilistic 
linear  model 


(L2)  X     =UY     _     +    (1  -  U   )A      /Wj,x 

v  n         n  n-S  n     n-(N+l) 

n 


for     n  =  1, 2, . . .  ,  where 


(1.3)  A     =  V  A     _    +   (l-V   )Y   . 

n         n  n-1       v       n     n 


The  model  of  (1.2)  and  (1.3)  will  he  termed  DARMA(l,N+l),  (discrete 
mixed  autoregressive-moving  average  process  with  autoregression  of 
order  1  and  moving  average  of  order  N+l) . 

If  we  start  the  process  with  A  /„ .  .  v  having  the  distribution 
tt  independent  of  {Y  ;  n  >  -N} ,  {U  } ,  {V  J,  and  {S  ),  then  the  X  's, 
n  i=  1,2,...,  will  be  shown  to  form  a  stationary  sequence  of  dependent 
discrete  random  variables  having  marginal  distribution  tt.  This 
stationary  sequence  is  in  general  not  Markovian,  although  it  will  be 
so  if    p  =  0.   Its  correlation  structure  is  determined  by  the  parameters 
p  and  (3  and  the  distribution  F,  Note  that  tt  can  be  any  distribution. 
Some  cases  of  discrete  distributions  of  particular  interest  are  obtained 
by  choosing  ir  to  be  geometric  or  Poisson. 

Certain  special  cases  of  the  MRMA(l,  N+l)  process  are  of 
particular  interest  and  their  consideration  will  make  the  nomenclature 
clear. 

(i)   The  MR(l)  process. 
If  P  =  0,  then 

(An-(N+1)-1    with  Probability  P> 
Yn-(N+1)      With  Probability  (1"P)« 


(A  ]  is  called  the  DAR(l)  process  (discrete  autoregressive  process  of 
order  1) . 


(ii)  The  DMA(N)  process 

If  3=1,  then  X  =  Y  _  =  Y  .   with  probability  F(k)   for 
'       n    n-S     n-k      *  J 

n 

k  =0,1, ...,N  where  F   is  the  distribution  of   S  .   In  this  case, 

n 

{X  )   is  called  a  DMA(N)  process  (discrete  moving  average  process  of 
order  N).   Note  that  if  {X  }   is  a  DMA(l)  process,  then 


n 


X  =  < 

n 


Y       with  probability  F(0), 


Y     ,  with  probability  1  -  F(0). 


(iii)     The  DARMA(l,  1)   process 

Finally,  if  N  =  0,  then  {X  }  will  be  termed  a  DARMA(l,l) 
process  (discrete  mixed  autoregressive-moving  average  process  both  of 
order  1)  with  parameters  p  and  p;  that  is, 


X  = 
n 


Y      with  probability  p, 


A  .    with  probability  (l-p) 


Note  that  the  DMA(l)  process  is  a  special  case  of  the  DARMA(l,l) 
process  when  p  =  0. 

(iv)   Independent  process 

If  p  =  1,  N  =  0  or  p  =  0,  p  =  0,  then   {Xj   is  a  sequence 
of  independent  random  variables  with  common  distribution  IT. 

The  model  of  (1.2)  and  (1.3)  is  really  the  backward  DABMA  model, 
The  forward  model  is  defined  in  a  similar  fashion.   However,  the  two, 
while  similar,  are  not  necessarily  equivalent.   This  is  because   (Xn) 


is  not  in  general  time  reversible  in  the  sense  that  {X..,...,X,  )  will 
not  in  general  have  the  same  distribution  as  {X  .,  X    ,  ...  ,  X  ). 

— it     —  K X  J.  —  1 

The  properties  of  one  model  can  be  derived  by  the  same  techniques  as 

those  of  the  other,  so  we  will  only  consider  the  backward  model. 

Note  that  the  MRMA  process  may  be  defined  using  any  sequence 

of  independent  identically  distributed  random  variables  {Y  } ,  not 

necessarily  discrete.  However,  (1.2)  and  (1.3)  show  that  {X  }  is 

n 

obtained  as  a  mixture  of  the  (Y  )   sequence.  As  a  result,  even  if 
the  distribution  of  Y   is  continuous,  a  realization  of  the  sequence 
{X  }  will  in  general  contain  many  runs  of  a  single  value.  This  seems 
to  be  the  major  drawback  to  using  this  scheme  to  obtain  a  sequence  of 
dependent  random  variables  with  a  specified  continuous  marginal  distri- 
bution and  correlation  structure.  Other  schemes  for  obtaining  sequences 
of  dependent  exponential  and  gamma  random  variables  have  been  proposed 
which  look  more  promising;  cf.  Lawrance  and  Lewis  (1977)*  Gaver  and 
Lewis  (1977),  and  Jacobs  and  Lewis  (1977). 

One  motivation  behind  the  DARMA  models  was  to  provide  a  simple 
scheme  for  obtaining  models  with  which  to  analyse  stationary  sequences 
of  dependent  discrete  random  variables  with  specified  marginal  distri- 
bution and  correlation  structure.   In  general,  there  is  not  much  beyond 
a  Markov  chain  model  which  is  overparametrized  for  statistical  purposes 
for  modelling  dependent  sequences  of  random  variables.  In  addition  it 
is  very  often  simple  to  show  from  data  that  the  correlation  structure 
of  the  sequence  is  not  Markovian.  The  DARMA  model  can  be  used  to  model 
nonMarkovian  sequences  of  discrete  random  variables;  an  observed 
sequence  of  this  kind  is  discussed  in  the  last  section. 

k 


Another  motivation  for  the  development  of  this  process  was  to 
provide  models  for  point  processes  in  which  the  data  is  given  in  terms 
of  counts  in  fixed  time  intervals  rather  than  the  exact  times  of  arrivals, 
Most  models  for  point  processes  beyond  the  Poisson  process  are  most 
easily  described  in  terms  of  times  of  arrivals  or  times  between  arrivals 
and  it  is  often  hard  to  obtain  results  concerning  the  joint  distribution 
of  counts  in  different  fixed  time  intervals.  We  feel  that  the  DAKMA 
models  will  be  of  use  in  such  situations.  There  is  also  the  possibility 
of  modelling  directly  the  binary  counting  process  in  a  discrete  time 
point  process.  This  is  discussed  in  Section  6. 


2,   SOME  PRELIMINARY  PROPERTIES  OF  THE  DARMA(l,N+l)  PROCESS 

In  this  section  we  will  give  some  properties  for  the  DARMA(l,N+l! 

process.  Unless  otherwise  indicated  we  will  assume  throughout  the  paper 

that  A  /N+,x  has  a  distribution  tt  and  is  independent  of  {Y  ),  {U  }, 

{V  ),  and  {S  }. 
1  n;        nJ 


2.1,  The  marginal  distribution  of  X  . 

We  will  first  show  that  X    as  defined  by  (1.2)  has  distribution 

tt  for  all  n.  To  this  end  we  note  from  the  expression  (1.3)  that  the 

random  variable  A   can  be  expanded  backwards  to  the  initial  value 
n 


A  /„  v  to  give  A  =  Y  .  with  probability  pJ(l-p)  for  0  <  j  <  N+n 
and  A  =  A  ,„ . ,     with  probability  p     ;  that  is,  A   is  a  mixture 
of  Y  ,  Y  i >  • • •  >  Y_N,  and  A  (N+1\ .  Hence,  for  i  in  the  state 

space  E 

N+n  . 

P(An=i}  =  Z  P(Vj=i}  pJ(l-p)  +  P{A-(N+l)=i}p  n   =  7r(i) 

for  n  =  -N,  -N+l,  ....   Similarly, 

N  N 

P{Yn   =i}  =  S  P{Y   =i,  S  =j}  =  ir(i)  Z     P(Sn=j}  =  ir(i) 
n"Sn     j=0    n"J     n  J=0    n 

for  i  €  E  and  n  =  1,2,...  .  From  (1.2)  and  (1.3)  it  now  follows  that 


P{Xn=i)  .  p  HY^^i)  ♦  (1-0)  P(An.(N+1)=l)  =  Ai) 


for  i  £   E  and  n  =  1,2,  ...  .  Hence,  the  marginal  distribution  of 

the  X  '  s  like  those  of  the  Y  '  s  is  it. 
n  n 


2.2.  Correlational  properties  of  {X  } 

Although  the  X  '  s  have  a  stationary  distribution  it,   the  X  '  s 

are  not  independent* as  are  the  Y  's.  This  can  be  seen  by  the  following 

calculation  of  the  covariance  between  X   and  X   ..  After  some 

n       n+j 

simplification 


(2.1)    E[Xn+JXn]  -  E[Xn+J]E[Xn] 

■  e2(E[Yn+d-S  +.  Vs  ]  -  E[Vi-S  +  3  ElYn-Sn]) 

0     n+0  n  °     n+j  n 

+  P(l-P)    I«I  An.(N+1)]    -  E[Yn+j.s        ]   E[An.(N+1)]) 

n+j  n-f-j 


+  P(l-P)    (E[An+..(N+1)   Yn.sJ    -  E[An+..(N+1)]   E[Yn.sJ) 

♦    (LP)2   (E[An+..(N+1)  An.(N+1)]    -  E[An+..(N+1)]   E[An.(K+1)]) 


The  covariance  of     Y     _       and    Y    ..  is 

n-Sn  n+1-sn+l 

=  A  PtSn+l*)  P(Sn=k+1)  Var  Yn-S 
k=0  n 


In  the  case  p  =  1,  (the  DMA(N)  model),  only  the  first  term  in  (2.1) 
is  nonzero  and  we  can  get  the  correlation  from  the  above  result. 
Putting  F(k)  =  P{S  =k)  we  have  for  the  correlation 

(2.2)    p^l)  =  =orr(Y    ,  I       ) 

n       n+1 


=  {E[Y  J._    _         Y     _   ]   -  E[Y  ±.    _        ]  E[Y     _   ])/Var  Y     a 
n+l-S    .  _      n-S  n+l-S   .-  n-S       '  n-S 

n+1  n  n+1  n  n 

N-l 
=    Z    F(k)  F(k+1)    . 
k=0 


By  similar  reasoning,  for  j  <  N 

N-j 
(2.3)     p^(d)  =  corr(Y    ,Y       )  =  Z  F(k)  F(k+j) 

n    J  n+j    k=0 

and  for  j  >  N,  PvtaQ)  =  0.  Note  that  these  expressions  do  not  depend  on 

n  and  thus  the  DMA(N)  process  is  second  order  covariance  stationary. 

We  will  now  compute  the  covariance  of  A   and  A  ..  which 
v  n       n-l 

appears  in  (2.1);  this  will  incidentally  give  us  the  correlation  structure 
of  the  DAR(l)  process. 

E[Vn-l]    "  E[V  E[Vl] 

"  ^VlVl]   "  E[Vl]  E[Vl" 

+   (1-p)   {E[Y  A     .]   -  E[Y  ]  E[A     ,]} 
v       '   l       n  n-lJ  nJ         n-lJ 

=  p  Var  A     , , 
n-l 


since  the  second  term  is  zero  because  Y   and  A  ,   are  independent. 

n       n-1         e 

This  is  because  A  n   is  a  function  only  of  Y  , ,  Y  „,....      By  an 

n-1  J      n-1   n-2' 

induction  argument  we  obtain  for  the  correlation  of  A   and  A  , . 

n       n+j 


(2.U)  PA(j)  =  corr(An,An+j)  =  pJ 


for  j  >  1.  Because  of  the  assumption  that  A  /„  ..  n  has  distribution  tt, 

(2.U)  does  not  depend  on  n  and  thus  the  autoregressive  process  is 

second-order  covariance  stationary. 

To  complete  the  result  for  the  general  DARMA(l,N+l)  process, 

we  compute  the  cross  covariances  between  the  sequences  {Y  _  }   and 

n— b 
n 

{A  } .  We  obtain 
1  nJ 


<2-5>      E[Yn+j-Sn+J  V(N+1)]  "  Et,Wsn+j}  EtV(^l)l  ■  ° 


for  i  >  1  since  S   .  only  takes  on  the  values  {0,1, ...,N}.  For 
°   —  n+j 

0  <  j  <  N 


<2-6'         EiYn-S    An-N+jl   "  E[Yn-S J  E!AnW 
n  n 

■  J  .  EfYn-kAn-N^   Sn=k]   "  *l\-rSrf*]  E[AnW 

rC(l-p)   F(N-j)    +   p(l-p)   F(N-j+l)    +   •••    +   pJ'(l-p)   F(NJ]   Var(Y1) 


For  j  >  N 


<2^     *Yn-Sn  V*J]  "  E[Yn-Sn]  EtAn-^] 


N 

=  Z  E[Y  .  A^  _.  .j  S  =k]  -  E[Y  .j  S  =k]  E[A     1 
.  -,    n-k  n-N+j'  n       n-k'  n      n-!H-jJ 


k=0 


=  [pJ'"N(l-p)  F(0)  +  pJ"N+1(l-p)  F(X)  +  ...  +  pj(l-p)  F(N)]var(Y1). 


Putting  everything  together  in  (2.1)  we  obtain  the  correlation 


of  X   and  X^. 
n       n+j 


°»M   "  (E[XnW   "  E[Xn]  «W,A*  Xl  • 


We  have 


p  N_1  p 

(2.8)  pjl)    =  ^     Z     F(k)   F(k+1)   +  P(l-P)  F(N)    (1-p)   +   (l-p)p 

k=0 


For     1  <  j  <  N 


(2.9)  (^(j) 


P  N"J  i-l_ 

=  P       Z    F(k)  F(k+j)  +  p(l-p)(l-p){F(N-j+l)  +  pF(N-j+2)  +  ...+  pJ  T'(N)} 

k=0 

+   (l-p)V    . 


10 


For  j  >  N+l 


(2.10)    p^j)  =  PJ"^  {0(l-p)(l-p)[F(O)+pF(l)  +  ...+pNF(N)]  +  (l-p)2  p1*1} 

Note  that  0  <  p^j)  <  1  and  for  j  >  N,  p^(j)   decreases 
geometrically  if  p  >  0  and  p  <  1.   Since  PwCd)   is  independent  of  n, 
the  MRMA(l,N+l)  process  is  second  order  covariance  stationary. 


2,3.  Invariance  under  transformations 

From  its  definition  we  note  that  X   is  a  mixture  of  the  random 

n 

variables  Y  ,  Y  ,,  ...  ,   Y  „  and  A_fN+1\;  i.e.,  it  is  a  random 

selection  of  one  and  only  one  of  these  random  variables.  Thus,  if  we 

transform  each  of  the  random  variables  Y  ,  Y  ,,...,  Y  „,  A  /„,, N  by 

n'   n-1'     '   -N'   -(N+l)   J 

the  same  function,  each  X   will  be  transformed  in  the  same  way  and 

its  distribution  will  be  that  of  the  transformed  Y  's. 

n 

Similar  remarks  apply  if  we  transform  the  X  's.   Note  that  in 
applying  a  common  transformation  individually  to  the  X  ' s  we  do  not 
affect  the  selection  procedure  and  therefore  the  correlation  structure 
of  the  transformed  process  is  the  same  as  that  of  the  untransformed 
process.  This  (marginal)  transformation  invariance  is  important  for 
statistical  analysis  of  the  process. 
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3.   THE  AUTOKEGRESSIVE  PROCESS  DAR(l) 

In  this  section  we  will  give  some  properties  of  the  MR(l) 
process  {A  } .  As  usual  we  will  assume  that  A  /   v  has  distribution  ir. 
By  the  results  of  Section  2,    {A  )   is  a  stationary  sequence  of  random 
variables  with  marginal  distribution  tt  and  correlations 


(3.1)  pA(j)  -  corr(An,An+j)  =  pJ  ,      J  >  1. 


The  spectrum  of  the  process  is  thus 


(3.2)  f(«D)  -s=JHv2  2  p.  (J)  cOS(o>j)}  =  ;±- 


2T     '  J-l  A  a,r  1  +  p2  -  2P  cos  * 


It  follows  from  (1.3)  that  {A  }  is  a  Markov  chain;  that  is, 
P{A    =  i|A  ,  ...,A  }  =  P{A    =  i|A  }  for  any  i  in  the  state  space  E. 
Further,  it  is  not  hard  to  show  that  the  transition  matrix  P  is 
given  by 

(3.3)     P{A    =  i|A  =  k}  =  p(k,i) 

p+(l-p)Tr(i)  for  k  =  i. 


Note  that  we  have  started  from  the  opposite  direction  from  that 
usually  taken  in  Markov  chain  theory;  we  have  specified  the  stationary 
distribution  associated  with  the  chain  first  and  specified  the  (Markovian) 
dependency  structure  by  a  single  parameter  p.  Moreover  changing  p 
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does  not  affect  it.      When  p  =  0  we  have  a  stationary  sequence  of 
independent  identically  distributed  random  variables  with  distribution  tt. 
The  fact  that  {A  }   is  a  Markov  chain  with  a  particularly  simple 
transition  function  P  makes  many  calculations  quite  easy.  For  example, 
in  discrete  time  series,  runs  of  given  values  of  the  random  variables 
X   are  useful  in  statistical  analyses.  Properties  of  these  runs  are 
easy  to  obtain  for  the  DAR(l)  process.  Thus  fix  a  state  i  £  E  and 
let  T.  =  inf{n  >  1:A  ^   i]  -  1;  Ti  is  the  length  of  a  run  of  i's 
starting  at  time  1,  where  length  can  be  0,1,...  .   Then 


Pf^  >  n}  =  P{A1  =  A2  =  ...  =  An  =  i}  =  ti  (i)  P(i,i)n": 


for  n  >  1  and  P{T.=0)  =  1  -  7r(i).  Thus 


(3.u)  «'i]- i-IIU-d-P^u)]  • 

If  p  =  0,  then  A  =Y    for   n  >  1  and  E[T,]  =  Tr(i)/[l-7r(i)  ] 
n    n         —  i 

as  expected  since  {A  }   is  a  sequence  of  independent  random  variables 
in  this  case.  Note  that  for  0  <  p  <  1 


e[t.]  >-  ¥(iL  ; 

L  i  -  1  -  tt(i) 


that  is,  the  expected  length  of  a  run  of  i's  for  a  DAR(l)  process  is 
always  greater  than  or  equal  to  the  expected  run  length  for  a  sequence 
of  independent  random  variables.  Moreover  the  inflation  in  the  expected 
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length  of  runs  is  uniform  for  all  states.  This  is  a  consequence  of 
the  fact  that  we  are  dealing  with  a  one  parameter  Markov  chain. 

It  is  also  not  hard  to  calculate  the  generating  function  for 
T . .  We  have  for  0  <  z  <  1 

*(z)  =  E  z^fT  =n}  =  [l-Tr(i)]  +  E  zn[P(i,i)  J11"1  [1-P(i,i)] 
n=0  n=l 

[l-Tr(i)][l-zp] 
"  1  -  Z7r(i)  -  zp[l-7r(i)]   ' 

Again,  if  p  =  0,  $(z)  reduces  to  -j   (jV  the  exP1"688*011  for  the 
generating  function  of  a  length  of  run  of  i  for  a  sequence  of 
independent  random  variables  with  marginal  distribution  ir. 
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THE  DMA  PROCESS 


In  this  section  we  will  consider  the  DMA.(n)  process  X  =  Y  „  . 

n    n-S 

n 

Note  that,  unlike  the  DAR(l)  process,  the  DMA(n)  process  is  not  MarkovLan 
in  general. 


U.l.   Correlation  properties. 

By  results  in  section  2,  {X  }  is  a  stationary  sequence  of  random 
variables  with  marginal  distribution  It     and  correlations 


N-j  N 

(l+.l)         PMA(J)  =  COrr(Xn,XnH)  =  ^  F(k)  F(k+J')  =  ^  F(v)  Y^~i] 

k=0  v=j 

for  1  <  j  <  N.  Also  PmU)  =0  for  j  >  N  and  Pm(o)  =  1.   Note 

that  when  N  =  1,  the  maximum  value  of  the  first  order  serial  correlation, 

max  p.,.  (l)  =  max  {F(o)  [l-F(o)]}=  1/k.  In  fact  one  can  show  that 

F(0)  m  F(0) 

for  any  N  >  1  the  maximum  first  order  serial  correlation  pWA (l) 
—  MA 

that  can  be  achieved  is  lA.   One  can  also  maximize  the  correlation 
at  any  point,  say  j,  by  making  F(j)  =  F(o)  =  1/2.   However,  all  the 
other  correlations  are  zero. 

For  the  spectrum  of  the  DMA(n)  process  we  have 
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oo 

(k.2)  f(o>)   =  ±-  [1  +  2     £     pMA(j)   cos(coj); 


27T  l  .    n     MA.  °  j 


then  if  we  define     p      (-j)   =  p     (j),  we  have 

MA       MA. 


(U.3)     fM  =  ^   £   elQ)J  £  ,F(v)  F(v-j) 


v=lo 


•  =  [(*  e"iCDj*  F(j))(E  eia*F(j))+l-  L(j)2] 
j=o         o=o  j=0 


N 
=  £=  [%(">)  <Pq(^)  +  1  -  £  F(j)2] 

27r  J-0 

=  ^[|cps(03)|2  *  1-  I  F(j)2] 

3=0 


where  cp_  is  the  characteristic  function  of  the  distribution  F  of 
the  random  variable  S.  Thus  we  can  model  a  broad  class  of  spectra 
f(oo).   if  F(o)  =  1  we  have  an  independent  identically  distributed 
sequence  and  a  flat  (constant)  spectrum. 

By  way  of  example,  it  is  worth  noting  that  we  have  restricted 
S  to  have  finite  support.   Then  (h.3)   is  a  polynomial  in  cos  oo  just 
like  any  moving  average  process.  The  finite  support  was  necessary 
to  allow  inclusion  of  the  autoregressive  tail  (1.2).   If  one  does  not 
want  to  add  this  tail,  then  there  is  no  reason  to  restrict  the  range 
of  S.   One  then  gets  a  much  broader  class  of  models  for  which  (U.3) 
in  particular  holds,  although  the  model  is  still  a  random  index  model. 
This  extended  model  is  not  as  broad  as  the  DARMA(l,N+-l)  model  in 
the  sense  that  one  cannot,  as  in  linear  models  (see  for  example, 
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Anderson,  1970 )  represent  the  tail  (1.3)  as  a  random  index  model  in 
which  the  random  indices  for  each  n  are  independent  random  variables, 

To  continue  with  the  example,  in  the  extended  moving  average 
model  let  S  have  a  geometric  distribution 


Then 


P{S=j}  =  F(j)  =  p(l-p)J        j  =  0,1, 


<p  (co)  =  P    , 

S      ,    ,,   \  ito 


1  -  (l-p)e" 


and 


27T 


P P .  -     y   2r,   ^2j 

[1  -  <l-p)e*>]   [1  -  (l-p).-*]   *  •  £  P  (1"P) 


27T 


5 +    gpd-p) 

[1+  (1-p)2  -  2(l-p)  cos  cd]   [1-  (l-p)2] 


The  initial  point  in  the  spectrum  is  related  to  the  amount  of 
long  term  dependence  there  is  in  the  process.  One  could  measure  this 
by  an  index  of  dispersion  (Cox  and  Lewis,  1966,  p.  71 ) 


Var(X  +  •••  +  X  ) 


k{E[X]}' 


and 


lim  .1  =  27T^L|  f(0+; 

k^co  K  E[xr 


For  the  moving  average  process,  f  (0+)  =  r—  [2  -  ^!  Q  ?(*■)    ] 
which  takes  values  between  1/27T  and  ir.      To  compare  the  moving 
average  process  to  the  DAR(l)  we  note  that  from  (3-2)  for  the  DAR(l^  process 
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f(00  =  ^.  [  (1-p2)/(1-p)2]  =  ^:  "(Yrjjy  which  is  always  greater  than  1/2tt 
if  p  >  0  and  increases  with  p  to  infinity.  Note  that  f (0+) 
for  a  sequence  of  independent  random  variables  is  1/2TT.  Thus  both 
the  moving  average  process  and  the  DAR  process  give  more  long  term 
dependence  than  a  sequence  of  independent  identically  distributed 
random  variables.   The  DAR(l)  process  allows  more  long  term  dependence 
than  the  moving  average  process. 

k.2.     Joint  distributions  and  time  reversibility 

Unless  otherwise  indicated  we  will  restrict  our  attention  to 
the  DMA(l)  process  in  the  remainder  of  this  section}  that  is,  if 
a   =  P{S  =0},  then 


(k.Q)  X 

n 


Y        with  probability  a 


Y        with  probability  (l-a)» 


In  this  case  PWA(l)  =  a(l-a)  and  P„A  ( j )  =  0  for  j  >  2. 

MA  MA  — 

It  is  not  hard  to  calculate  the  joint  Laplace-Stieltjes 

transforms  of  the  joint  distributions  of  random  variables  in  the 

-sX 
DMA.(l)  sequence  but  it  is  tedious.  For  example,  if  r(s)  =  E[e    ], 

then  from  (k.Q) 

*2(sl,S2)  =  Ef^Pf-SiV  S2X2}] 

o 

=  (l-a)  r(s  )  r(s  )  +  ad-a)  rCs^s  )  +  a  t(b  )  r(s  ) 
=  r(sx)  r(s2)  [l  -  a(i-a)]  +  a(i-a)  r(s1  +  sg). 

18 


Similarly,   by  conditioning  arguments  we  obtain 

t3(si,S2'S3)   -  E[exp  f-s^   -  s2X2   -   s^  }] 

=  (l-a)  r(s1)  ^(s^s^)  +  a(i-a)  r(s1  +  s  )  r(s  ) 

+  a  (l-a)  r(s1)  r(s2  +  s  )  +  a5  r(s1)  r(s  )  r(s  ) 
=  r(s1)  r(s2)  r(s  )  [1  -  2ad-2)] 

+  [r(s1)  r(s2  +  s3)  +  r(sL  +  s2)  r(s3)]  a/(i-a) 


and 


VWS3'S^   =  E[exP{-siXi  "  S2X2   "  S3X3  '  SA}] 

=  (l-a)  r(sx)  ^(s1,s2,s^)  +  a(i-a)  r(s1  +  sg)  *2(s,,  s^) 

+  a  (l-a  )  r(s1)  r  U2+s5).  r(s^)   •  cr  (l-a)  rU^  r(s2)  r(s  ,+s^) 

+■  ar(s1)  r(s2)  r(s,)  r(s^) 

=  r(s]_)  r(s2)  r(s3)  r(s^)[i-3a(i-a)  +  a2(i-a)  -a5  (l-a)] 

+  r(s1)  r(s2)  r(s3+s^)[a(i-a)  -a2(i-a)  +  a3(i-a)] 

+  r(s1)  r(s2  +  s  )  r(s^)  a(i-a) 

+  r(s1  +  s2)  rCs^)  r(s^)  [a(l-a)  -  a2(l-a)  +  ex5 (l-a)] 

+  r(s1  +  s  )  y(s    +■  s^)  a  (l-a) 

One  interest  in  the  joint  distributions  of  the  random  variables 
is  to  look  at  the  time  reversibility  of  the  process.   One  reason  for 
concern  with  time-reversibility  is  the  following.   The  EMA(l)  process 
(exponential  moving  average  of  order  l)  of  Lawrance  and  Lewis  (l977) 
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is  not  time  reversible  even  though  this  fact  cannot  be  determined 

from  second-order  properties  of  the  process.   Consequently  one 

cannot  distinguish  between  a  and  (l-a)  in  the  spectrum  of  the 

EMA1  process.   However,  by  using  higher  order  moments  it  is  possible 

to  distinguish  between  a  and  (l-a).   For  the  DMA(l)  process  the 

fact  that  p  .(l)  =  a(l-a)  means  we  cannot  use  it  to  distinguish 

between  a  and  (l-a).   The  time  reversibility  for  the  DMA.(l)  process 

would  mean  that  we  might  not  be  able  to  distinguish  between  a  and 

(l-a)  even  by  using  higher  order  moments. 

Since  \|r  (s_,s0,...,s  )  =  t  (s  ,s  ,,...,8.)  for  n  =  2,3A> 
n  12      n     n  n  n-1      1 

it  seems  likely  that  the  DMA.(l)  process  is  time  reversible.   In  order 

to  show  time  reversibility  we  need  to  show  that  t  (s, ,s_, . . . ,s  )  = 

J  n  12      n 

ty  (s  ,s  ...... s.)  for  any  nonnegative  s,,...,  s  .   For  simplicity 

n  n  n-1      1  1      n         r  * 

consider  the  terms  a  =  r(s,  )  r(s~  +  s,)  r(si,  )  ...  r(s  )  and 

12    3    4        n 

b  =  r(s1)  r(s2)  ...  r(sn_5)  r(sn_2  +  s^)  r(sn)  of  *n(i1,«2,...,.n). 

In  the  expression  for  t  (s, ,s_,...,s  )  the  term  a  has  a  coefficient 
r  n  1  2      n 

of 
n 
£  PfS^O  or  1,  S2=0,  Sy   +1,  S^=  1,...,S  *1,  S.+1-0,...,Sn=0] 

J  —  j 

n 

=  E  pfs^o,  s2=o,  s3=  +i,  sk=  +i,...,s  =  +i,  sJ+1=o,...,sn=o} 

n 
+  Z  P{SL=  +1,  s2=o,  S^=  +1,  S^=  +1,...,S .=  +1,  s   =o,...,sn=o} 

since  the  event  associated  with  the  term  a  is  that  X_  and  X,  pick 
the  same  Y.  and  that  all  the  other  X.'s  pick  distinct  Y.'s  from 
each  other.   Similarly,  the  term  b  has  a  coefficient  of 
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n-2 


Z  P{Sn=0  or  1,  S   -1,  Sn_2=0,  S   =0,...,S  =0,  S   =1,...,S  =1] 

i=l  J 


Since  the  S  's  are  independent  and  only  take  the  values  0  and  1 

and  the  Y  '  s  are  independent,  the  coefficients  of  the  terms  a  and  b 

of  ty  (s,,srt,...,  s  )  are  equal. 
n  1  2      n 

Similar  arguments  can  be  used  to  show  that  the  coefficients 
of  the  terms 


h'1  W1 

n    r(s.)  r(s.  +s        )       n     r(s      .)  r(s    ^s        ) 

i=l  x  Jl     Jl  i=2  Jl  J2     J2 


^k'^k-l"1  n 

n       r(s.        . )  r(s    +  s        )  r(s  ) 

i=2  Jk-1  Jk       Jk         i=j   +2 


and 


Jx-2  j2-jr2 

n       y(s     .        .)  y(s     ,  ..  +  s        )      n       rU     .    . 

i=i        ^^-r1      n-Jk+1     n-Jk    i-i        n-Jk-x 


of  \|r  (s,,...,s  )   are  equal  for  any  sequences  of  integers 
n  1      n 

l  <  dx  <  j-l+i  <  d2  <  J2+i  <  o3  <  •  *  •  <  dk  <  Jk+1  <  n 

and  k  >  1.   Thus  ^(s^.  •  •  ,sr)  =  ^n(sn>- • • >  s]_)   for  all  n. 
Hence  the  DMA(l)  process  is  time  reversible. 
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k.3.      Run  lengths  for  the  DMA(l)   process. 

We  will  now  consider  length  of  runs  for  a  DMA(l)  process. 
Fix  a  state  i   in     E     and  let     T.    =  inf{n  >  1:X  7^i}-l     the  length 
of  a  run  of     i     initiated  at  time  1  where  length  can  be  0,1,...    . 

We  will  first  compute     E[T.].      Let     a^  =  1     and     a     =  P{X  =X  =   •• 

*  i  0  n  L   1     2 

for     n  >  1.      Then 


=  X  =1) 

n 


&1  =  PfX^i) 


.2  _/,  n2 


a2  =  PfX^x  =i}  =   (l-a)  Tr(±)a1  +  a(l-a)  7r(i)a0  +  a    7r(i) 


and  hy  induction 


Vl  =  P(Xl=i'    -    '   Xn+l=i) 


=   (l-a)  ir(iU    +  a(l-a)  Tr(i)a     ,  +  a  (l-a)  ir(i)  a    _  + 
n  n-l  n-2 

+  a  (l-a)  ttu;  &0  +  a        ttu) 


Thus 


E[T.]   =     E     a 

1  i      n 

n=l 


n=l 


00  00 

=  ir(i)  +   (l-a)  ir(i)     E    a     +  a(l-a)  ir(i)  :       ■;,  \       E    an 

n=l  1-ottUJ         , 

+  °    T{i)2  T^TiT  +  a(l-a)  7r(i)  i^FTiT  ao 


=  ir(i)  + 


a 


ir(i)     +  cflKi)  -  a  ir(i) 


l  -  air(i) 


(l-a)  ir(i)  + 


ad-a)ir(i) 


l-arr(i) 


E[T.  ] 


Solving  the  last  equation  for  E[T. ]  we  obtain 
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(h   Q^  V\T    1    -      P(i).     _  tt(j)[I   +   Q(l-u)(l    -   7T(lj)1 

{h-9)  E[Ti]   -  l^flT  -  [1  -rr(i)][l  ■  a(l-a)  ttUH 

where 

(^.10)  p(i)   =  Tr(i)[l  +  a(l-a)]   -  7r(i)2  a(l-ctf 

=  7r(i)   +  a(l-a)  -rr(i)[l  -  7r(i)]    . 

If  a  is  either  0  or  1,  then  {X  }   is  a  sequence  of  independent 
random  variables  and  E[T.]  =  7r(i)/(l-Tr(i) )  as  expected.  Note  that 
E[T.  ]  >  Tr(i)/[l-7i(i)]  for  0  <  a  <   1;  that  is,  the  expected  length 
of  a  run  of  i  for  a  DMA(l)  process  is  greater  than  the  expected 
length  of  a  run  for  a  sequence  of  independent  random  variables. 
For  a  given  distribution  tt,  the  maximum  value  for  E[T.]  occurs  when 
a  =  l/2.  In  this  case 

We  now  turn  our  attention  to  the  generating  function  of  T. . 
Fix  j  ^   i  in  the  state  space  E  and  let 

b  -P{X  -  •••  -Xn-1,  Xn+1  =  J),   n>l 


and 


bQ  =  P{X1=j}  . 


Using  an  induction  argument  we  obtain  for  n  >   1 
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bn  =  (1-a)  Tr(i)  b^  +  an-(i)    (l-a)^  +  o?7r(i)2(l-a)bn_3 

+  ...  +  an-177(i)n-1(i-a)b    +  cftr(i)n  cen-(j) 


Thus 


J    z\  =  b    +    £    A  fz(l-a)7r(i)  +  (l-a)z   f    anTr(i)nzn} 
n=0  n=0         n  n=l 


00 

n  ,1  Nn  n 


+  Qf7r(j)    £    a  7r(i)    z 
n=l 


0       n^         n)  v  1  -  cwrUTz     [  1  -  az7r(l 


After  some  simplification  we  obtain 


(lull)  •(z)   =    t     ACT^n)   =       Mi)l[l-  z^)  g&SU 

n=0  1-Z7r(i)  -  z  tt(1)  [  l-r;(i)  ]  a(l-a) 


Note  that  for  a  =  0  or  1,  $(z)=[l-ir(i)]/[l-Z7r(i)  ]  as  expected. 
Higher  order  moments  of  the  run  lengths  can  be  obtained  from  (4.11) 
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5.   THE  BINARY  DARMA(l,l)  PROC 

In  this  section  we  will  consider  a  DARMA(l, 1)  process  in 
which  X   takes  only  the  values  0  and  1;  that  is, 


' 


(5.1)  X  =< 

\S  !  n 


Y      with  probability  p, 


<.-■ 


and 


A  =( 

n   \ 


A      with  probability  (l-p) 


A      with  probability  p, 


Y      with  probability  (l-p), 


where  {Y  }   is  a  sequence  of  independent  random  variables  taking  the 
values  0  and  1  with  common  distribution  it.     Note  that  the  DARMA(l,l) 
process  is  not  Markovian  in  general. 

Time  series  of  binary  random  variables  are  of  particular 
importance  for  modelling  the  differential  counting  process  in 
discrete  time  point  processes.  Klotz  (1973)  and  Kanter  (1975) 
have  given  a  model  which  is  different  from  the  binary  DARMA(l,i^ 
process. 

Setting  N  =  0  in  (2. 8) -(2. 10)  we  obtain  the  correlation 

of  X   and  X   . ,  j  >  1, 
n       n+j    — 


_  J'1 


p^  =   corr(Xn,Xn+J)  =  p^pd-p)  (l-.V  +  {!-, 
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The  spectrum  of  the  process  is  thus 


(5.2) 


f(o>)  =  ~  {1  +  2  E  qJL3)    cos(aij)} 
3=1 


_  1  |2c((3,p)  [co 
"  24        1  " 


s  (a  -   p]  -  2p  cos  op  +  p' 


where 


2p  cos  to  +  p 


c(p,p)  =  p(l-p)(l-p)  +  (l-$)2p  . 


^ 


We  will  now  consider  some  properties  of  lengths  of  runs. 
For  fixed  i  €  {0,1},  let  T.  =  inf{n  >  1;X  ^i)  -  1  as  before. 
We  will  calculate  E[T.]  and  the  generating  function  for  T..  To 
begin,  note  that  although  {X  }  is  not  a  Markov  chain,  {(A  ,X  )j  n=l, 2, . 
is  a  Markov  chain.  For  i,  I,    j,  k  in  {0,1) 


.} 


P{An+]=d,  Xn+1  =  k|An=i,  Xn=^}  =  0^1,3 


is 


independent  of  I.     Letting  Q,   denote  the  matrix  whose  (i,j)  entry 
0^1,3)  we  have 

p(l-p)  +  [l-p(l-p)]  tt(0)      (1-P)(1-P)tt(1) 

P(i-p)  tt(o)  Pp  ir(o) 


Q  = 


and 


Q  = 


3p  ir(l)  p(X-p)  tt(1) 

(l-p)(l-P)  tt(0)    (l-p)p  +  [l-(l-p)p]  ir(l) 
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Note  that  P{Tn  >  n|A.  =  i)  =  q"(I,0)  +  Q?(i,l)  s 


QS(i'E: 


for  i  =0,1.   Hence 


E[T0|AQ=i]  =  £  Qj(i,E)  =  RQ(i,E)  -  1 
n=l 


where  1^(1,  j'1  =£  _o  Qn(i>j)   with  1=0^  heing  the  identity  matrix 
and  R  (i,E)  =  R  (i,0)  +  R  (i,l).   It  is  not  hard  to  show  that 


R0  =  (I  "  Q0 


-1 


where 


(i-Q0) 


-l   l 


"l  -  Pp  tt(o) 

p(l-p)  77(0) 


(1-3)  (l-p)  77(1) 
1  -  p(l-p)  -  [1  -  p(l-p)]  tt(0) 


with 


A  -  det(l-QQ)  =  [1-tt(0)]{1-p(1-3)[1-M0)]-M0)[1-P(1-p)]}  • 


Thus 


(5.3: 


E[TQ] 


tt(0)  RQ(0,E)  +  tt(1)  Rq(1,E)  -  1 

7r(0)fl  -  gp  +  e(l-p-frf2gp)[l  -  tt(0)]} 


[l-TT(0)]fl-p(l-3)[l-p  tt(0)]  -  pV(0J[l-P(l-p)]} 
Tr(0){l+corr(Xn,Xn+1)+pp-p)  +  tt(0)2  (-corr(Xn,Xn+1) 


-  2Pp  +  p] 


X-p(l-P)  +  7r(0){-l+2p-3eP-corr(Xn,Xn+1)}  +  7r(or{corr(Xn,Xn+1)  +  23p-d; 

7r(0){corr(Xn,Xn+1)+l-p+pp-7r(0Hcorr(Xn,Xn+1)  +  2pp-p)] 
[1  -  7r(0)]{l-p+pP-7T(0)(corr(Xn,Xn+1)  +  2{3p-p)} 


after  some  simplification. 

Similarly  one  can  show  that 
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(5.U) 

ELT1J  "  [l-77(lj  ](1-P(l-P)[l-Ml5]  -  Mi)[i-P(i-P)]} 


7T(1)C1  +  corr(Xn,Xn+1)+pP-p}  +  71  (l)2{-corr(Xn,Xn+1)-2(3p  +  p} 
l-p(l-P)  +7r(l){-l  +  2p-3pp-corr(Xn,Xn+1)}+7](l)2{corr(Xn,Xn+1)+2Pp-p} 

tt(1)  {corr(Xn,Xn+1)+l-p+pP-7r(l)(corr(Xn,Xn+1)  +  2Pp-p)} 
[l-7r(l)]{l-p+£p-77(l)(corr(Xn,Xn+1)  +  20p-p/]  ' 


Note  that  the  expected  length  of  a  run  of  i  for  the  binary  MBMA(l,  1) 

process  is  always  greater  than  or  equal  to  the  expected  length  of  a 

run  of  i  for  the  independent  case. 

We  now  turn  our  attention  to  the  computation  of 
00 

$-U)  =E  _^  z^PfT.sn}  for  i  =  0,1  and  0  <  z  <  1.  To  begin,  note  that 

1  *-'n=0  1  —       — 

P{T0=n|A0=i)    =ZQ°(i,j)[l  -  C^(J,E)]    . 

J 
Thus 

00  00 

£    z^CT  =n|A =i}   =  ££     zWiJ)    [l-Q0(j,E)J    . 
n=0  j  n=0 

It  is  not  hard  to  show  that     Y,°°  =$  ^^n   =    ^  "   ZSd^~   '  wnere 


,  .         fl  -  zpp  tt(0)  z(l-p)(l-p)   tt(1)  "I 

(l-zo   )"1  =  -7L- 

A{z)        zp(l-p)  ir(0)  1  -  z{p(l-P)  +  [l-p(l-P)]  7r(0)}J 


with 
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A(z)    -  1    -   zp(l-p)    +   ztt(O)    {-p-l+p(l-p)) 

+    z27,(0)[l-7r(0)]{-a(l-p)    +  2pp(l-p)}    +   z2n(0)2ep   . 

After  some  manipulation  we  obtain 
•0(«) 

=  l-zp(l-p)+T7(0)[l-7;(0)]{zfi2-z32p-zp-zp2p)    +  7?(0)[zp-l]  -zgp  TTJQf 

l-zp(l-p)+Z7r(0){-pp-l+p(l-p)}+z277(0)[l-7r(0)]{-p(l-p)+2p3(l-p)]+z2 

In  a  similar  manner  one  can  show  that 


&x(z) 

=     l-zp(l-p)  +  7r(l)[l-7T(l)](zp2-zp2p-zp-zp2p)  +  tt(1)[zp-1]  +  zpprr(l)2 
l-zp(l-p)-tZ7T(l){-pp-l+p(l-p))+z27r(l)[l-7T(l)I-p(l-p)+2pp(l-p))+z27T(l)2ep 


Higher  moments  for  the  runs  can  be  obtained  from  the  generating 
functions.   These  are  important  in  determining  to  what  degree  the  binary 
DARMA(l,l)  process  differs  from  Markov  models,  e.g.  the  DAR(l)  model, 
the  differentiation  considered  here  being  to  what  extent  the  distribution 
of  run  lengths  departs  from  a  geometric  distribution. 
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6.   TELEPHONE  ERROR  DATA 

We  discuss  here  very  briefly  a  case  in  which  the  binary  DARMA(l,l) 
model  may  be  of  use;  in  particular,  we  do  this  to  illustrate  some  of  the 
formulas  given  in  the  previous  section.  Another  recent  example  of  the 
need  for  discrete  time-series  models  is  Gaver,  Lavenberg  and  Price  (1976) 

Cox  and  Lewis  (1966,  p.  173)  discussed  data  consisting  of  errors 

which  occurred  during  the  transmission  of  binary  data  over  a  telephone 

line.   Let  X  =  1  indicate  that  the  nth  transmitted  bit  is  in  error, 
n 

while  X  =  0  indicates  that  no  error  occurred.   Models  which  postulate 
that  bit  errors  occur  independently  or  according  to  a  Markov  chain  such 
as  the  binary  DAR(l)  process  predict  that  the  runs  of  ones  and  runs  of 
zero  will  both  have  geometric  distributions.   But  the  runs  of  zeros  are 
just  the  intervals  (or  number  of  bits)  between  errors  without  the  times 
consisting  of  1  bit  between  errors,  and  these  were  shown  by  Berger  and 
Mandelbrot  (1963)  and  Lewis  and  Cox  (1966)  to  be  nongeometric.   In  fact 
they  are  highly  skewed,   long-tailed  distributions  which  led  Berger  and 
Mandelbrot  (1963)  to  postulate  a  model  in  which  intervals  between  errors 
were  assumed  to  be  independent  with  Pareto-type  distributions. 

The  problem  with  the  Berger-Mandelbrot  model  was  that  the 
intervals  between  errors  were  found  to  be  dependent  (Lewis  and  Cox,  1966) 
Moreover  there  are  a  disproportionate  number  of  1-bit  between  error 
intervals;  128  out  of  673  intervals,  while  the  longest  interval  between 
errors  is  85 > 993  bits.   This  suggests  that  modelling  the  binary  X 
process  may  be  a  better  approach  than  modelling  the  intervals,  although 
the  modelling  process  must  be  nonMarkovian. 
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The  binary  DARMA(l,l)  process  is  a  candidate  process  for 
modelling  this  process;  in  particular  one  would  like  to  know  whether 
the  runs  of  zeros  for  p  between  zero  (Markovian)  and  one  (independent; 
produce  highly-skewed  run-length  distributions.   The  question  is  too  broad 
to  be  considered  here,  involving  also  estimation  of  suitable  values 
of  (3  and  p,  and  will  be  considered  elsewhere.   Here  we  will  examine 
only  the  effect  of  p  on  E(T  )   and  the  plausibility  of  the  model. 

For  the  error  data,  672  bits  out  of  1,106,148  transmitted  were 
in  error,  so  that  we  can  estimate  tt(1)  as 

T?(l)  =  l  -  tt(0)  =  x  xo^S  =  0.0006075  • 

Thus  from  (3.4)  with   p  =  0  we  compute  that  the  expected  lengths  of 
runs  of  zeros  and  ones,  given  that  they  occur,  are  respectively 
1/tt(1)  =  1,645.09  and  l/rr(0)  =  I.OOO608;  the  observed  values  from  the 
data  are  1,911.27  and  1.235,  both  much  longer  than  predicted  under 
independence  assumption  (p  =0). 


(5.3). 


In  Table  1  we  give  values  of  E(T  )  computed  from  the  formula 


Note  that  for  small  p,  e.g.  p  =  0.1,  the  value  of  E(TQ) 


first  increases  with  increasing  p,  and  then  decreases.  This  is  character- 
istic behavior  for  the  process  when  it  is  almost  a  moving  average.  For 
large  p,  we  find  E(T0)   decreasing  with  p.   In  particular  p  has  a 
large  effect  on  E(T  );  it  remains  to  be  seen  how  p  effects  the  whole 
distribution  of  runs. 
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The  estimated  first  and  second  serial  correlation  coefficients 

for  the  data,    p(l)      and     p(2)      are  0.1^0  and  0.121  respectively.      This 

i 

is  consistent  with  the  models  restriction  to  positive  serial  correlation 

coefficients.   From  the  expressions  for  PM(l)   and.  PM(2)   given  in 
Section  5  we  see  that  p(2)/p(l)  -0.6k     is  a  rough  estimate  for 
thus  p  is  relatively  large  for  this  data.   With  the  proviso  that 
is  relatively  large  it  might  he  possible  to  find  unique  values  of  p 
and  p  which,  with  the  estimated  Tr(l),  would  make  (5-3  '  and  (^.k) 
equal  to  the  estimated  E(T  )  and  E(T  ).   (This  is  not  possible  in 
general  since  the  expressions  (5-3)  and  (5 A)  are  not  single-valued 
functions  of  (3  for  small  fixed  p.  )  An  alternative  is  to  use  the 
estimate  of  p  and  E(T  ),  with  ir(o)  and  p(l),  in  15-3)  and  solve 
for  p.   The  rough  estimate  obtained  this  way  is  p  -  0. 

It  appears  to  be  possible  to  estimate  (3  and  p  in  a  more 
systematic  way  using  higher  order  joint  moments  of  the  X  's.  This 
will  be  discussed  elsewhere. 
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